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I. INTRODUCTION 


Since the first nonlinear generalization of the Dirac equation by Ivanenko [T], the nonlinear Dirac (NLD) equation 
has emerged naturally as a practical model in many physical systems, such as extended particles mg, the gap 
solitons in nonlinear optics [g , light solitons in waveguide arrays and experimental realization of an optical analog for 
relativistic quantum mechanics IMS], Bose-Einstein condensates in honeycomb optical lattices [g, phenomenological 
models of quantum chromodynamics nni, as well as matter influencing the evolution of the Universe in cosmology 
mi. To maintain the Lorentz invariance of the NLD equation, the self-interaction Lagrangian can be built up from 
the bilinear covariants. Different self-interactions give rise to different NLD equations. Several interesting models 
have been proposed and investigated based on the scalar bilinear covariant HMS], the vector bilinear covariant US], 
the axial vector bilinear covariant m, both scalar and pseudoscalar bilinear covariants m, both scalar and vector 
bilinear covariants [HI |20| among others. 

A key feature of these NLD equations is that they allow solitary wave solutions or particle-like solutions - localized 
solutions with finite energy and charge [H]. That is, the particles appear as intense localized regions of field which 
can be recognized as the basic ingredient in the description of extended objects in quantum field theory [22j . For the 
NLD equation in (1-1-1) dimensions {i.e. one time dimension plus one space dimension), several analytical solitary 
wave solutions were derived for the cubic nonlinearity [23l [2g , for fractional nonlinearity [Sg as well as for general 
nonlinearity uniiiiiiTi by using explicitly the constraints resulting from energy-momentum conservation; and this 
is well summarized by Mathieu [5g. With the help of the analytical expressions of the NLD solitary wave solutions, 
the interaction dynamics among them has been studied and rich nonlinear phenomena have been revealed in a series 
of works [29H8g . 

An interesting topic for the NLD equation solitary waves is the stability issue, which has been the central topic in 
works spread out over several decades that is still ongoing. Analytical studies of the NLD solitary wave stability face 
serious obstacles [55H57] , while results of computer simulations are contradictory nzuMsg. Numerical results inferred 
that both the multi-hump profile and high-order nonlinearity could undermine the stability during the scattering of 
the NLD solitary waves [301 |3g. For the NLD equation with scalar-scalar interactions (i.e. the Soler model) the 
solitary wave solutions can have either one hump or two humps. Quite recently, for the Soler model, we found that 
all stable NLD solitary waves have a one-hump profile, but not all one-hump waves are stable, while all waves with 
two humps are unstable |41j . This result is in agreement with the rigorous analysis in the nonrelativistic limit [42j . In 
order to further understand the behavior and the stability of NLD solitary waves, the NLD equation in the presence 
of external potentials has been investigated |43H46j and a sufficient dynamical condition for instability to arise was 
postulated through a collective coordinates (CC) theory [jg. In this work, we will continue to study the NLD solitary 
waves under external forces. 

For the forced nonlinear Schrodinger (NLS) equation when subject to an external force of the form f{x) = 
r exp{—iKx), the authors found [37H43] that intrinsic soliton oscillations are excited, i.e., the soliton amplitude, 
width, phase, momentum, and velocity all oscillate with the same frequency. This behavior was predicted by a collec¬ 
tive coordinates theory and was confirmed by numerical simulations. Moreover, one specific plane wave phonon (short 
for a linear excitation) with wavenumber k = —K is also excited such that the total momentum in a transformed NLD 
equation is conserved. This phonon mode was not included in the CC theory and had to be calculated separately m- 

In the present paper we consider the relativistic generalization of our previous work on the forced NLS equation, 
namely the behavior of solitary wave solutions to the NLD equation when subjected to an external force which is now 
a two-component spinor. In Sec. |ll|we review exact analytical solutions for the unperturbed NLD equation. In Sec. 
Ill we present the NLD equation with external force fj{x,t) = rj exp[i(iyjt — Kjx)], j = 1,2, and the corresponding 
Lagrangian density. Using the energy-momentum tensor we show that the total energy is conserved if the force is 
time independent {vj = 0). 

For the case Ki = K 2 = K, Vj — 0 and zero dissipation we perform in Sec. IV a transformation such that the 
transformed NLD equation is invariant under space translations and thus the momentum is conserved. In Sec. jV] 
we make a variational ansatz with three collective coordinates. All integrals that appear in the Lagrangian can be 
performed exactly and we finally have a set of 3 ODEs as CC equations. For a special case these equations can be 
simplified and an approximate analytical solution can be obtained. Solutions are also obtained in the non-relativistic 
regime, when AT = 0, by an expansion up to order where v is the velocity. 

In Sec. VI the spectrum of the linear excitations (phonons) is calculated and together with the numerical solutions 


of the CC equations this is compared with the results from our numerical simulations (Sec. VII). We always obtain 
periodic solutions and the spectra of these solutions exhibit two dominant peaks: The phonon frequency VIk = 
and the intrinsic oscillation frequency flsim- The frequency flsim agrees nearly perfectly with the prediction 
rice from solving numerically the CC equations (Tablejl]). The soliton position q(t) performs small oscillations around a 
mean trajectory Vsim t- This translational motion is only weakly affected by the intrinsic soliton oscillations. Further 
Vec agrees with Vsim within an error of about 14% (Table |T|. The reason is that the plane wave phonons with k = —K 
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are not taken into account in the CC theory. The summary of our main results is contained in Sec. VIII 


II. REVIEW OF EXACT SOLUTIONS TO THE NLD EQUATION 


In this section we review the exact solitary wave solutions to the NLD equation, 


- m)^'+ = 0 , ( 2 . 1 ) 

where we use the representation for the 1+1 dimensional Dirac Gamma matrices: 7 ® = 173 ; 7 ^ = ia 2 , which we also 

used in |46| . The solitary wave solution in the rest frame is represented by 

Vh(a:, t) = e-“‘^(x) = e—‘ ) , (2.2) 

where A and B satisfy 

HA 

+ {m + uj)B- g^{A^ - B'^YB = 0 , 
ax 

^ + {m-uj)A-gYA^ - B^'A = 0 . 
ax 

(2.3) 


The solutions of these equations vanishing at infinity are 


A = 

B = 


1 (m + Lo) cosh^(s:/3x) 

(k + l)/3^ 

m + u} cosh(2K/3x) 

_g"^{m + u) cosh(2K,dx))_ 

(« + l)/3^ 

/ (m — w) sinh^(K/3x) 

m + u} cosh(2K/3x) 

_g'^{m + UJ cosh(2K/3x))_ 


(2.4) 


where /3 = . We are interested in bound state solutions that correspond to positive frequency w > 0 and 

which have energies in the rest frame less than the mass parameter m, i.e. uj < m. 

Because of Lorentz invariance we can find the solution in a frame moving with velocity v with respect to the rest 
frame. The Lorentz boost is given in terms of the rapidity variable 7 as follows (here c = 1): 

1 V 

n = tanh 77 ; 7 = —== = coshr;; sinh? 7 =—==. (2.5) 

V 1 — V 1 — 

In the moving frame, the transformation law for spinors tells us that: 


cosh(r 7 / 2 ) sinh(??/ 2 ) \ f «')’[ 7 (x - nt), 7 (t - ux)] \ 

^ ^ ^ sinh( 7 y/ 2 ) cosh(? 7/2 J ^'^[ 7 (x — nt), 7 (t — ux)] y ’ 


( 2 . 6 ) 


since 


cosh(?7/2) = V(1 + 7)/2; sinh(77/2) = 1/(7 - l)/2. (2.7) 

This in component form reads: 

'I'i(x,t) = (cosh(7y/2)A(x') + isinh(77/2)i?(x')) , 

'I' 2 (x,t) = (sinh(77/2)A(x') + icosh(? 7 / 2 )i?(x')) , (2.8) 

where 


x' = 7 (x — vt); t' = 7 (t — vx). 


Note that cosh^(? 7 / 2 ) + sinh^(? 7 / 2 ) = cosh 77 = 7 . 


(2.9) 
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III. EXTERNALLY DRIVEN NLD EQUATION 

In previous papers [111 HI] we investigated the externally driven NLS equation 


+ SjIj = re 


(3.1) 


where /x is the dissipation coefficient, and r and K are constants. This equation can be derived by means of a 
generalization of the Euler-Lagrange equation 


d dL d dC 


dC 


dF 


dt dipt dx dipt dip* dipt 


where the Lagrangian density reads 


^ - 'tptP’) - \'>Px? + + 5\ip\^ - ^ 

I n I 

and the dissipation function density is given by 

F =-i^i{iptip* -ipti’)- 

For the NLD case we instead consider a two-component spinor forcing term 


^iP* - re'^^ip 


^ V hlx,t) ) 


with the NLD equation 

(*7^9^ - m)4' -h = f{x, t) - . 

In what follows we will generalize our choice for the NLS equation by choosing 

j = l,2, 


(3.2) 


(3.3) 


(3.4) 


(3.5) 


(3.6) 


(3.7) 


with real parameters Xj, Vj and Kj. Note that the phase of / is invariant under Lorentz transformations. As the 
second component of the spinor 'I' is the so-called “small component”, which is smaller than the first component by 
the factor a = \/ (m — uj)/{m + ix), we will only consider cases with r 2 = e?”!, where e = 0{a) or smaller. 

Equation (3.6) can be derived in a standard fashion from the Lagrangian density 


^ = J - f-^ + Co{b), 

. ^ ) “T 1 


(3.8) 


where Co{b) is determined later on and b = lima;_>±oo 'L(a;,t). The term in the Lagrangian density which pertains to 
forcing can be written as 


£3 = -2i?e(/>), 

and the full interaction part of the Lagrangian density is now 

Cl = - T/ - />. 

tv ~r -L 

The generalized Euler-Lagrange equation can be written as 

dC _ dF 

^'"did^^) ~d^~ d{dt^ ’ 
where the dissipation function density is now 

F = - 5ti'7°4'). 


(3.9) 


(3.10) 


(3.11) 


(3.12) 
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The adjoint equation comes from the Euler-Lagrange equation: 


dC. dC dT 

from this we get the adjoint driven NLD equation 

— = /+ • 


(3.13) 

(3.14) 


To generalize our discussion of external forces from the NTS equation to the NLD equation we have included a 
dissipation term in our general formulations. However, in most sections that follow we will concentrate on the case 
where the dissipation coefficient /x = 0, so that the energy is conserved. 

When Ki = K 2 = K, Vi = 0 and ^ = 0 we choose KL = irn (n is an integer, 2L is the total length of the system) 
and periodic boundary conditions such that fj{x,t) goes to rj and the wave function goes to a constant spinor b for 
L —>■ 00 . In that case a constant solution of the forced NLD equation with 4' = 6 satisfies the nonlinear equation: 

— mb + g‘^(bb)'^b = r, (3.15) 


where r represents a spinor with components ri and r 2 . 


A. Energy flow equations and the conservation of energy 


From the NLD equation with external sources and the definition of the energy-momentum tensor: 




we have that 


where 


= F'' 


The energy density is given by 


where now 


F’' = ^(aV) + (5V)«'. 


TOO = - _ [4-719,,4' - 9„4-714-] -f m4-4- - £/ - £0, 


£, = ^(§4-)"+1-/4/-§/, 

Av “r J- 


(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 


and £0 is chosen so that vanishes at a: = ±£, when L —>■ 00 . Therefore, from Eqs. (3.15) and (3.19) we obtain 

(3.21) 


£0 = mbb - ^—-{bb)'^~^^ + br + rb = —mbb + - - — -{bb)'^~^^, 


g‘^{2K F 1 ) - 


J1 

K -|- 1 ' ^ ' ' ' K -|- 1 

Now we will assume that in the lab frame f{x, t) is independent of t and of the form: 


fj{x) = rje j = l,2. 


with real parameters rj and Kj. In that case from Eq. (3.18), we have that = 0 and 


dtT^^ -b 9„Tii’ = 0, 


where 


r™ = — [4-7193,4- - 93,4-714-] -I- m4-4-- - —(4-4-)'^+! + 4-/ -I- /4- 

2 Kj —)— 1 


7^'° = "2 ■ 


(3.22) 

(3.23) 

(3.24) 

(3.25) 


Integrating Eq. (3.23), and under the assumption that T^°{+oo,t) — T^^{—oo,t) = 0, then the energy of the driven 
NLD equation, 

r+oo 

dxT°°, (3.26) 


^total _ 


is conserved. 
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IV. TRANSFORMED NLD EQUATION WITH EXTERNAL FORCE 


Let us consider the case of external force (3.5) and (3.7) with Vj = 0 and Ki = K 2 = K, 


fix) 


ri 

r2 




(4.1) 


After the following transformation 


Eq. (3.8) becomes 




C = 



[Xl^d^,x - ^^,X1^X] - ixi-xx + -xr-rx + Kxi^X + ^o(a), 

A€ J. 


(4.2) 


(4.3) 


where the constant Co{a) now depends on a constant vector a = lima;_j.±oo x(2;, t), determined by the algebraic 
equations: 


— ma + g^{aa)'^a = r — K^^a. 


(4.4) 


From the Euler-Lagrange equation: 


^^d{d^x) dx 

we obtain the following perturbed NLD equation 

(*7^5^ - m)x + ff^(xx)''x = x- K-f'^x- 

The corresponding adjoint equation reads 

- idf.xi^ -mx + ff^(xx)'"x = x- Kxi^. 


(4.5) 


(4.6) 


(4.7) 


Equation (4.6) is not only invariant under time-translation, but also under space-translation, 
energy and the total momentum should be conserved quantities. Indeed, multiplying Eq. 


Therefore, the total 
(4.6[) from the left by Xt 


and Eq. (4.7) from the right by yt, and then adding both expressions, we again obtain the continuity equation (3.23), 
now with 


2 

= “b [xi^^^X- dxXl^x] +mxx -^(XX)'"'^^ + Xr + rx - Kxx^x - ^oia), 

^ z K H- i 

where £ 0 ( 0 ) is chosen as 


(4.8) 


jCqM = rnaa -(aa)"^^^ -|- ar -|- fa — Kaj^a, 

K + 1 


(4.9) 


so that T™ vanishes at a; = ±L, when L ^ 00 . Moreover, 


T^° = -'o[xa^x-xi\t]. 


(4.10) 


Integ rating Eq. (3.23), and again assuming T^^{+oo,t) — T^°(—oo,t) = 0, the energy of the driven NLD equation 
(4.4) given by 




r +00 


dxT^^, 


(4.11) 


is conserved. Inserting the transformation (4.2) into the equation (3.26), it can be verified that 


rptotal Tpiotal 

— h/^ 


(4.12) 
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we obtain the continuity equation for momentum flow 


Now multiplying Eq. (4.6) from the left by Xx and Eq. (4.7) from the right by Xx^ and then adding both expressions, 

(4.13) 


_yOi _|-2^11 ^ Q 

dt ^ dx ^ 


with 


and 


= 7 ) [Xxl^X - Xl°Xx\ , 


'^x = d [x7°Xt - Xt7°x] - mxx + Kxi^X + -^ixxT""^ -xr-rx + ^o(a). 

^ Z Av H- i 


(4.14) 


(4.15) 


Integrating Eq. (4.13), if T^^(+oo,t) — r^^(—oo,t) = 0, then the momentum of the transformed driven NLD equation 
is given by 


r+oo 


^x = 


dx- [xxl°X - Xl^Xx] , 


which is also conserved. Now inserting the transformation (4.2) into (4.16) 

2+00 

Px = 


dx \ - ^'7°^x] } ■ 


(4.16) 


(4.17) 


The right hand side of Eq. ( |4.17 ) cannot be separated into two integrals due to the fact that 4'^^ 4' does not vanish at 
X —)■ ± 00 . However, in simulations when we deal with a finite domain, we have 


where now 



— KQ + 

(4.18) 


p-\-L 



Px = 

L 

dx^ [Xxl°X- Xl°Xx] , 

(4.19) 


p-\-L 



Q = 

L 


(4.20) 


p-\-L 



p^ = 

L 

dx- [ 4 ',, 7 ° 4 ' - 4 ' 7 ° 4 '^] . 

(4.21) 


V. VARIATIONAL (COLLECTIVE COORDINATE) ANSATZ FOR THE NLD EQUATION WITH 

EXTERNAL DRIVING FORCES 

Our ansatz for the trial variational wave function is to assume that because of the smallness of the perturbation 
the main modification to our exact solutions to the NLD equation (without driving forces) is that the parameters 
describing the position q{t), width parameter j3{t) and phase become time dependent functions. We assume 
that the driving term is specihed in the lab frame, and that the initial condition on the solitary wave is that it is a 
Lorentz boosted exact solution moving with velocity v. To describe the position of the solitary wave we introduce the 
parameter q{t) which replaces vt for the unforced case. We then let the width parameter /3 and thus u) = rn? — /3^ 
become functions of time. We next rewrite the phase of the exact solution as 

ujt' = 7 a;(t — vx) —)■ (j){t) — p{t){x — q{t)) (5-1) 

to mimic our parametrization of the collective coordinates in the nonlinear Schrddinger equation. Next, we let 
p(t) = L>j{t )x{q )q be determined from u}{t) and q{t) and let the phase 4>(t) be an independent collective variable. That 
is, in Eq. (|2.8|) we replace 


vt —>■ g(t); (3 —>■ /3(t); ivt' = jujit — vx) —>■ (f>(t) — p{t){x — q{t)), 


(5.2) 












where p{t) = 'y{t)oj{t)q{t). 

Thus our trial wave function in component form is given by: 


^-1 (cc, t) = (cosh + i sinh , 




<{x,t) = (sinh |a(z) + i cosh 


(5.3) 


where z = coshry (a: — q{t)). Note that w, which was a parameter in Eq. (2.4), now is time dependent because 


of w = — /3^(t). Using the trial wave function Eq. (5.3) we can determine the effective Lagrangian for the 

variational parameters. Writing the Lagrangian density as 


C — Cl + C2 + 


where 


£1= 1 , 


£2= —mdttl/ 
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K + 1 




K+1 . 


C3 = -^f-f'if. 


Integrating over x and changing integration variables to 0 one obtains 

Li = j dxCi = Q (^pq + ^ — p tanh — Iq cosh p — Jq tanh p, 

where the charge 

Q = J dz[A'^{z) + B'^{z)], 

and the rest frame kinetic energy Iq = Hi 

Iq = j dz {B'A - A'B) = Hi, 


are given by Eqs. (Al) and (A2), respectively, in the Appendix. Here B'{x') = ^ , and 


dx' ’ 


and 


(5.4) 


(5.5) 


(5.6) 


(5.7) 


(5.8) 


= J dz (ha - AH^ , 

(5.9) 

dA dA dz 

^ _ :__ :__ 

dt dz dt ’ 

[x — q{t)) cosh? 7 , we have 

(5.10) 


— = —q cosh p — z tanh pp 


(5.11) 


and 


Jq = — cosh pqlo — p tanh p J dz^z {AB' — BA'). 

The integrand in the second term is odd in z, so the integral vanishes and we are left with; 

Li = J dxCi = Q (pq + (j) — p tanh p^ — Iq (cosh p — q sinh p ), 


(5.12) 


(5.13) 


£2 = / dxCo = — 


m 


-I 1 + 


coshry (K+l)coshr 7 ’ 


(5.14) 
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where 


Ii= J dz{A^{z) - B‘^{z)) = J ~ 

For L 3 we have 

L 3 = -2 J dxRe [fi^i - / 2 ^' 2 ] = ~J ^^^3- 


K+1 hi -\- 1 K 1 

= - ^-^3 = —o --^ 1 - 


g^K 


(5.15) 


(5.16) 


We obtain for the integrand 


£3 = — 2 ri cos(^ — Kiq) I cosh ^A{z) cos 

I 2 7 


p + Ki 


P + Ki ■ u'd \ ■ 

-z — sinh —n[z) sin- 




+ 2 r 2 cos (0 — K 2 q) { sinh ^A(z) cos ^ z — cosh ^B{z) sin ^ z 

2 7 2 7 


(5.17) 


where we have not included terms that are odd in z. Performing the integration we get 

L 3 = —2— cos{(j) — Kiq) ^cosh ^ Ji — sinh — 2— cos(^ — K 2 q) ^sinh ^ J 2 — cosh ^-^ 2 ^ , (5.18) 


Jj{uj,q) = J dzA{z) cos 
Nj{oj,q)= / dzB{z)s\\i 


P + Kj 
1 


• cos bi 


z = 


g^/u! cosh a jir’ 


7 

P + Kj 

2/37 ’ 


a, = bi = a, cosh ^: 


i/u), 


p + Kj TT sin bj 


7 


g\/iij cosh GjiT 


(5.19) 

(5.20) 


The integrals Ii, I 2 , Jj and Nj are done exactly in the Appendix. Putting all terms together and using the fact that 
q = V = tanh 77 we obtain: 


L= <30 — Iq sechg — mli sechg + 


9 


K + 1 


I 2 sechp + L 3 , 


2 -k ( riCOs{(j) - Kiq) r 2 cos{cj) - K 2 q) 


^2 r, 


^ gj^/oj \ coshoiTT ^ cosho27r 

h} Tj rj 7 ] 

Cj = cosh - cos bj — sinh - sin bj , Sj = sinh - cos bj — cosh - sin bj, j = 1 , 2 . 


(5.21) 


Since we are using the exact solutions of the NLD equation as our trial wave functions for the forced problem, the 
integrals Iq, Ii and I 2 are related since for the NLD equation without the presence of external forces, the solitary wave 
with u = 0 obeys the relationship [23 




K + 1 


(^^)«+i = 0. 


(5.22) 


For our problem this converts into 


a;(A^ + B^) - m{A^ 


B"^) 


K + 1 


(A 2 


^ 2 )«+l ^ 


(5.23) 


Integrating this relationship we obtain: 


mil — 


I 2 — ujQ — H 2 - ojQ — 0. 


K + 1 


Using this relation to replace Ii and I 2 in L we obtain 


L = Q<j>- -{lo+ujQ) - U{q,q,P,<p) = - — - U{q,q, I3,(j)), 

7 7 


(5.24) 


(5.25) 


where U = —L 3 , and Mq = /q + u}Q is the rest frame energy of the solitary wave for k = 1. 
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From Eq. (3.12) we can calculate the dissipation function F for the CC equations. We find 

r+oo 


/ -\-oo 

dxlm{^>^ dt'^) 

-OO 


' —OO 

f*+oo 


= 2fi 


dz 


sinh ?7 {AB — BA) — cosh 77 {pq + + B^) 


cosh ?7 . 

We recognize the integrals as being related to Jo = — coshTygJo and Q, so we obtain 


(5.26) 


F = —2fi Iq sinh ijq + Q{pq + (j)) 


(5.27) 


We can simplify this by introducing the boosted rest frame mass: 


M = 7 M 0 = 7 (/o + ujQ) (5.28) 

and use the definition of p{t) = jujq so that 

F = -2pL{M<f + Q0). (5.29) 


This is the relativistic generalization of our expression that we found for the forced NTS equation 
ready to derive Lagrange’s equations for the collective coordinates using Eq. (5.25). From 


Now we are 


d dL dL _ dF 

dt dq dq dq ’ 

we obtain 

Jt 

where 

dt dq dq dq' 

We also have a contribution from dissipation from the equation 

d dL dL _ dF 

dt dij) 9(j) d<j )' 

which gives us a first order differential equation for uj 

dU 

Q = Q'{uj)uj = -2pQ - —, 


(5.30) 

(5.31) 

(5.32) 

(5.33) 

(5.34) 


where the prime denotes the derivative with respect to uj. 

As L does not depend on /3, the final Lagrange equation is dL/djd = 0. After changing to the variable w = i/m^ — /32 
we have 



(5.35) 


This leads to a first order differential equation for (jj 


Q'{uj)^ 


-M'iuj) 


dU 
duj' 


(5.36) 


In what follows we will make the simplification /r = 0 in our comparison of the CC equations with the numerical 
solution of the forced NLD equation, so that we have energy conservation as a check for our simulations. 
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A. Simplification at k = 1, Ki = K, ri = r, r 2 = jJ, — 0 
In that case U = —L 3 simplifies to be: 

U{q, q, OJ, 4>) = 2r cos{(j) — Kq) 




cosh ?7 


coshry 


where 


, TTCOsb , ,, 7rsin6 

•/i(w,g) = — -j= -r—; ^i(w,g) = 


g-y/w cosh ATT ’ 


5 y/aJcosh OTT 


and 


V^K oj^q + K I— - - 

2/37 2/37 

b{u},q) = acosh~^(m/a;), 7 ( 9 ) = coshr]{q) = (1 — 


From the expressions (Al) and (A3) in the Appendix we obtain 

Mq{uj) = ujQ'{u}). 


Inserting Eqs. (5.40), (Al), and (5.37) in Eqs. ( 5.36| ) and (5.34) we obtain, respectively, 

• OJ g^r(j(uj)oj"^ cosicj) — Kq) \coshrj/2 dJi{q,oj) sinh 7//2 97Vi(q, w) 

(j)= - 

7 




OJ = -- 


g^r(3oj^ sm{(j) — Kq) 




cosh rj doj cosh t] doj 


cosh rj 


cosh ?7 


In the special case of /r = 0, from Eqs. (5.31) and (5.32) we obtain 


d ,, ^ \ d dU dU 


(5.37) 

(5.38) 

(5.39) 

(5.40) 

(5.41) 

(5.42) 

(5.43) 


B. Solutions when K = vq = qo = 0 


When we look in the rest frame where r; = 0 and look for solutions, we notice when K = vq = 0 that a = b = 0. 
Thus A^i = 0 and 


TT / ! \ dJ ]_ I TT 


(5.44) 


Therefore, the equation of motion given by (5.43) is always satisfied. The Eqs. (5.41) and (5.42), respectively, become: 

gTTrfdoj^/'^ cos (j) 


= 00 + 


OJ = -- 


2m? 
'Kgrfioj^l'^ sin()) 




Note that li (jjq and we start with (jjq = mr (n is an integer) we have 

( 77 rr/ 3 a;^/^(— 1 )” cos ^ 


= U! + 


OJ = -- 


2m^ 

7 r( 7 r/ 3 a;^/^(— 1 )” sin 0 


This has an unstable stationary solution for n = 1, namely 


= 0 , ()) = 0 , Wo = 


g-KT 

2m^ 


/3ov^- 


For small girr 


Wo = 


2 2 2 
g TT^T- 

4m^ 


(5.45) 

(5.46) 

(5.47) 

(5.48) 

(5.49) 

(5.50) 


This solution is unstable to small perturbations. 
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0(t) -w(g)t 



FIG. 1. Left panel: — io{g)t vs t, Right panel: Lij{t) vs t. 

solving numerically the CC equations. Parameters are: cjq = 


ai(t) 



Blue curves are the analytic result. Red curves are the result of 
1.9, (j(0) = 0, A" = 0, ri = 0.01, r 2 = 0, m = 1 and g — 1. 


1. Analytic solutions for (j>o = 7r/2, K = v = 0 and g-Kr jrnf ^ 1 


If we choose (/)o = 7’‘/2, we get the simple equations: 


= w — 


U) = -- 


sin(^ 
2 m^ ’ 

TTgrjSuj^^'^ cos ^ 
m? 


One can expand these equations as a power series in g = girr/m?. To order g^ we obtain: 


m = I + 


Wo + y - 2wo) 


j ^ 3g/3„|cosy) - 11 ^ .. .inyt) ^ ^ ^ 


2\/wo 


8 wo 


u}{t) = Wo - gy^lio sin(wot) + g'^ sin^(wot) “ 2wo^ + O (g^) , 


(5.51) 

(5.52) 

(5.53) 

(5.54) 


where /3o = \/ rn? — lOq. We notice that the term in that is linear in t shows a constant shift from its initial value 
Wo, namely wo —t w(g) = wo + ^g"^ (m^ — 2wo). For g = m = 1, r = 0.01 and wo = 0.9 we find w(g) = 0.899694. For 
these initial data we have compared our analytic solutions with the numerical solutions of the CC equations. At all 
times the w(t) found analytically from Eq. ( |5.54| ) tracks the numerical solution almost perfectly slightly getting out 
of phase at late times. The phase </), after we eliminate the linear growth as determined numerically, gradually starts 
diverging from the analytical result Eq. (5.53) even at modest times. These results are shown in Fig. 


C. Variational method for the transformed NLD equation with external force, Ki — K 2 = K 


Let us now consider the case of an external force (3.7) with iq = 0 and Ki = K 2 = K. From (4.2) and the ansatz 


(5.3) it follows that 


where 


Xi 


{x,t) = (cosh^A(z) +Isinh^B(z)) 


2 

V 


2 

V 


(5.55) 


X 2 ix,t) = [sinh jA{z) +i cosh jB{z)j 


p = jujq + K, 
(j) = (j) — Kq. 


(5.56) 


Inserting (5.55) into the Lagrangian density (4.3) and integrating over x we obtain 


L = Q{^ + Kq)-^-Uiq,l3,^), 


7 


(5.57) 
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where 


U{q,P,^) 


2tt cos (j){riC — r2S) 
cosh ATT 


(5.58) 


with a = p/(2 /37), b = acosh”^ m/w, and C and 
not depend on q because the potential U in Eq. 
From the Lagrange equation 


S are defined by Eq. (5.21) with h instead of b. 
(5.58) no longer depends on q. 


Remarkably, L does 


d dL dL 
dt dY dY 


with Y = q, it is clear that the canonical momentum 


(5.59) 




— = jAIoq + KQ 




is conserved. Moreover, we obtain 


{^ + Kq) 


dQ 

duj 


1 dMo dJJ 

7 duj dio 

dQ dtj 
dt Q(p 


(5.60) 


(5.61) 

(5.62) 


from (5.59) with Y = ui, Y = (p, respectively. Multiplying Eq. (5.60) by q, using (5.61)-(5.62), and integrating over t 
we obtain that the energy 


^cc = — + (Pq - KQ)q + 17, 


(5.63) 


is conserved. 


Inserting the ansatz (5.55) into the expression for the field momentum (4.16) we obtain 


Px =PQ + qiHi. 


(5.64) 


q must be a conserved quantity. In simulations 


This means that in the solutions o f col lective coordinates equa tions. Pc 
of the transformed NLD equation (4.4), the momentum (4.16) must be conserved and equal to 7’x(0). 


A similar remark holds for the energy. Indeed, inserting the ansatz (5.55) into (4.11) and integrating over x we 
obtain 


Ex = 1^2 + U[q,l3,p). 


(5.65) 


This energy is eq ual to the energy of the driven NLD equation E^, obtained by inserting the ansatz (5.3) into ( |3.26| ). 
The energy Ecc (5.63) obtained from the solutions of the collective coordinates equations must be constant, while 
the energies E^*^ and computed from the simulations of the driven NLD equation and transformed NLD 

equation, respectively, not only must be equal to each other, but also they must be constant if the energy density 
satisfies certain condition at the boundaries. 

As U i n Eq. (5.58) does not depend on q{t); using v{t) = q as a. new collective coordinate we obtain one algebraic 
equation (5.60) and the following two differential equations of first order: 


sin^ 

'yrrP cosh ott 


{riC - r2S) ^ - 




~ UJ 

(j) = —Kq H- 

7 


dU 

2nP duj ’ 


(5.66) 

(5.67) 
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where in 7 we should also replace q = v{t) and 


duj 

QhjJ 

a 

dU 

dq 

dg 

C 4 

S 4 


u rr~ ^ ~ , 27rcos^(riC'^ - r 2 S'^) 

- TTuaui tanh an H-=-^-. 

2uj 57 ^^ cosh ttTT 


Kuj + jm^q 


27/33 


buj = du, cosh- 


_i m am 


oj tup, 


-b. 

-k 


TJ ~ Tl 

sinh - cos b + cosh - sin h 
2 2 

cosh — cos b + sinh — sin b 
2 2 . 


27r cos (j) 
g^y/Zj cosh an 
uj — Kq"f 
2/3 ’ 


(-( 97 ^ + ndqtai]hna){riC - r 2 S) + riCg - r 2 Sq) , 


^ 1 1 ^ 
- Qn cosh - 

Tl ~ V ^ 

sinh - cos b + cosh — sin b 

2 ^ w 

2 2 J 

1 m 

- an cosh — 

7? ~ 7? -■ 

cosh - cos b + sinh - sin b 

2 ® UJ 

L 2 2 J 


(5.68) 

(5.69) 

(5.70) 

(5.71) 

(5.72) 

(5.73) 

(5.74) 

(5.75) 


The Eqs. (5.66)-(5.67) can be obtained using the Lag range equat ion ( 5.59) and setting Y = r espec tively. Now 
we need to solve two first order differential equations (5.66) and (5.67 ) and one algebraic equation (5.60), where the 
three unknowns are v{t) and Then, using the substitution (5.56) one can obtain ip{t). 

In the above equations the subscripts on the variables d, b, C, S refer to (partial) derivatives with respect to the 
subscript variable. 


D. Non-relativistic (NR) regime with K\ = K, r\ = r, r 2 = 0. 


In the non-relativistic regime, u = g <C 1. Therefore, 

sinhSiT] ^ Siq; cosh —>■ 1 ; p ^ coq; 7 —>■ 1 , 

where in our case i5i = 1/2 or (5i = 1. Thus in determining U, 


d{q,u}) = = d cosh ^{m/uj). 


(5.76) 


(5.77) 


For determining the NR equations for p and oj we need an expression for U valid up to terms linear in q. We find 


U = 2r cos ( 


- 1^1 ( 9 , w) 


Thus the equations of motion (5.66) and (5.67) become, respectively, 

g‘^rl3uj‘^ sincf) 


(jj = -- 




-^ 1 ( 9 ,w) - -Ni{q,u) 


= —Kq + uj — 


g'^rPu'^ cosp 


dJi{q,Uj) _ q dNi{q,uj) 
duj 2 duj 


and the canonical momentum (5.60) becomes: 

2KI3 4m 


P — 

J g — 9 


tanh 


m — UJ . dU 
m + UJ ^ dq 


(5.78) 

(5.79) 

(5.80) 

(5.81) 


g^uj g- 

To determine the conserved Pq up to terms linear in u, we need an expression for U up to terms of order We find 

tzK 


17 = -- 


•cos(</)sech( ^^/;f_^, ) 
^g^/uJ {m^ — 


(oo -I- aiv a2V^ + ...), 


(5.82) 
















































where we have 
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.2 .,2^3/2_J ^cosh 1 (^) 


oq = —8 (m^ — cos 


y 2\/m? — 


U}^ 


(5.83) 


and 


W 2 2n 3/2 . /Kc0sh-1(^)\ ■ /^^COSh-l((^) 

ai = 4 sin' -1 i , ,2^ ^ I - 


/ A COSn I-t) \ . / ■) 2\ 1-1 ■ / A COE 

, + 4a; - o;^) cosh ^ — sin —^ 

y 2 Vto 2 - a ;2 y ^ \^> \2^t 


+ 47 ra; irr? — oj^) tanh ( - 


ttK \ /Xcosh-i(^) 


COS 


2\/w? — ui^ J I 2\fn^^^^u? 


Note that when K = 0, ai —>■ 0. We will display 02 for K = 0 below. 


(5.84) 


1. Nonrelativistic solutions with K = 0 


When it' = 0, it is easy to expand U up to order and obtain the two equations of motion and one constraint 
equation up to order v. One has: 

27rr cos{(j)) 


Uni - 


gV^ 


X i - 


— a;2 [3m? + (tt^ — 3) a;^) + — oP- cosh ^ + 2a; [m? — o;^) cosh ^ \ 

8 (w? — a; 2 )^^^ 


(5.85) 


We have 


dUnr 27rr cos{(j)) 
dv gy/uJ 

V — a;^ ( 3 to ^ + (tt^ — 3) o;^) + — o;^ cosh“^ {w? — o;^) cosh“^ (^)) 

4 (m^ — 

The non-relativistic conservation law for P„ when K = 0, valid to linear order in v is then 


Pq = Mo{u})v - 


dUn 


(5.86) 


(5.87) 


This allows us to explicitly solve for u as a function of a;(t), (j>{t) and the initial values of v,(j> and to. Note that 
because the velocity corrections to U start at v^, the leading terms to solve for (j) and uj are the same equations that 
we needed to solve in the 1 ; = 0 case, namely: 


UJ = -- 


= U! + 


ngrPuj^^'^ sin (j) 


gTTrPuj^/'^ cost/) 

2m^ ’ 


(5.88) 

(5.89) 


which we have already solved explicitly as a power series in g. These solutions are given in Eqs. (5.53), (5.54). We 
have that 

V 

— = PqivOi^oAo) 

Vo 


Mo{oj) 


2TTrcos{(j)) ( 3 ™^ + (tt^ - 3) a;2) + a;2/3(a;) cosh ^ 2a;/32(a;) cosh ^ (^)) 


gy/uj 


4 / 33 ( 0 ;) 


-1 


(5.90) 
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v(t) 



FIG. 2. w as a function of time t for an initial non-relativistic velocity vo = 1/50. Parameters are: cjo = 0.9, K — 0, 
ri = r = 0.01, r 2 = 0, m = 1, g — 1. The upper curve is the analytic solution. 


For the case (/)o = f, one has Pq{vo,uJo,(j)o) = Mq(ujq). For the initial conditions we considered earlier when vq = 0, 
namely K = 0,r = 1/100, m = l,g = 1, and initial velocity vq = 1/50 we get the result shown in Fig. In that 
figure we compare the analytic result with the solution of the CC equations. The analytic result is slightly higher at 
the maxima of v{t). 


VI. SPECTRUM OF THE LINEAR EXCITATIONS (PHONONS) 


Similar to the case of the forced NFS equation [^, the external force fj = [see Eq. (3.22)] excites not 

only soliton oscillations, but also a plane wave phonon (short for a linear excitation) such that the total momentum 

is conserved. _ 

The general solution of the linearized NLD equation without damping [Eq. (3.6) with /r = 0] and Ki = K 2 = K, 
reads 


with the phonon dispersion curve 




i^phik) = \/m2 + 


( 6 . 1 ) 


( 6 . 2 ) 


and arbitrary, but small b. In the case ri = r and r 2 = 0, which was considered in Sec. V, the spinor c has the simple 
form 


02 

“if 


—m 

K 


= \/ m? + 


(6.3) 


We choose \K\ m, then the second component of c is much smaller than the first one [see below Eq. (3.7)]. 

In the presence of a soliton we need to know the phonons only far away from the soliton. For example, when both 
the soliton and the phonon move to the right, at the far left the phonon wave function is given by Eq. (6.1) with 
0 = 0, while at the far right there is a phase shift 9^0. 

We now calculate the charge Qph and the momentum Pph of the phonon by integrating T'flE' and over the 
interval —L < x < L. Here we must distinguish two cases: k ^ —K and k = —K. In the latter case the x-dependent 
parts in the integrands of Pph and Qph drop out. We then obtain 


Qph = 2L{b^b + c^c) + 2L 
We also find Pph = —KQph so that 


(6.4) 


Pph + KQph — 0 . 


( 6 . 5 ) 
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The left hand side of Eq. ( |6.5[ ) is the total momentum in the transformed system discussed in Sec. IV. It is zero 
because the wave function describing the phonon, 


Xph = = b e 




C, 


( 6 . 6 ) 


is simply a homogeneous oscillation and does not travel in the transformed system. 

In our simulations which will be presented in the next section, we use a finite system of length 2L and periodic 
boundary conditions. The predicted phonon mode with wavenumber k = —K is clearly identified in the spectrum of 
Q{t), see panel (e) of Figs. and The observed frequencies W 2 agree very well with given by Eq. (6.3) for 


K = QimdK = with L = 100. 


The phonon mode is also seen indirectly in the spectrum of the maximum of the charge density p{x, t) shown in 
panel (f) of Figs. 00 and|^ This is a local quantity which is used for the computation of the soliton position q{t), 
in contrast to the global quantity Q(t) which is obtained by integration over the whole system. The phonon frequency 
a ;2 is observed in the difference W 3 = a ;2 — wi, where uji is identified as the frequency of the intrinsic soliton oscillations 
discussed in the next section. 


VII. SIMULATIONS VS NUMERICAL SOLUTIONS OF COLLECTIVE COORDINATES EQUATIONS 

In order to obtain the numerical solutions of the collective coordinates equations we need to set initial conditions for 
q(0), a;(0) and (/)(0). In simulations we use the soliton solution of the unperturbed NLD equation with the same initial 
conditions. We would like to stress here that arbitrary sets of these initial conditions produce different quantities for 
Pq and (see Fig. |^. Numerical solutions of collective coordinates must conserve Pq and Ecc, whereas in simulations 
and E-^ are conserved. Therefore, good agreement between simulations and numerical solutions is expected only 
for the initial conditions that guarantee Pq = P^ and Pec = E^. The simplest case is to choose initially 0(0) = ± 7 r/ 2 , 
g(0) = 0 and arbitrary values for q{0) and a;(0) (see Fig.|^. 






FIG. 3. Left panel: Pq (solid line) and (dashed line) vs initial phase. Right panel: Ecc (solid line) and (dashed line 
superimposed) vs initial phase. The maximum difference between E^c and is of order of 10~^. Parameters are: ujo = 0.9, 
( 7 ( 0 ) = 0.1, K = — Stt/IOO, ri = r = 0.01, r 2 — 0, m = 1, g = 1. 


The CC theory leads to the algebraic equation ( 5.60 ) and the two DDEs (5.66), and ( 5.67| ) which are solved by 
a MATHEMATICA program. The driven NLD Eq. (3.6) is a PDE for which various numerical schemes have been 
proposed that are reviewed in Ref. [SS]. It is also reported there that the operator splitting (OS) method performs 
better than the other schemes in terms of accuracy and efficiency. Therefore we have applied a fourth-order OS 
method in a recent paper on the stability of solitons in the undriven NLD equation m and the readers are referred to 
Refs. [33] and [H] for a detailed description of the method. For the driven NLD equation we again employ the same 
scheme used in Ref. [41], but instead of nonreflecting boundary conditions, we take periodic boundary conditions. 
This has the advantage that tests of the conservation laws for momentum and energy (see Secs. IV and III) yield 
an accuracy of the order of 10“^°and 10“^"^, respectively. Here we adopt the computational domain [—100,100], (i.e. 
L = 100), the time step At = 0.025, and the final time = 800. 
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FIG. 4. Left panel: Pq = vs initial velocity. Right panel: i5, 
K = —Stt/IOO, n = r = 0.01, r 2 = 0, m = 1, g — 1. 



= vs initial velocity. Parameters are: ljq = 0.9, (?i(0) = n/2, 


In the previous section we have identified plane wave phonons in the spectra of the soliton charge and amplitude 
(maximum of the charge density). Now we discuss the peaks due to the intrinsic oscillations of the soliton shape 
and velocity and compare with our CC results. The highest peak wi = flsim in the spectra of Q{t) is always close 
to the initial value uiq = 0.9 and always agrees nearly perfectly with the predicted frequency flee of the CC theory, 
see Figs. and fm and Table |Tj The dependency on the parameter K is weak which means that the periodicity 
of the force fi = rexp(—ziFx) has little influence on the intrinsic oscillations; this includes the case K = 0 where 
the force is homogeneous. Exactly the same frequencies are observed in the spectra of the soliton amplitude [in the 


simulations this is the maximum of the charge density, and in the CC theory this is a = 2 ^™ Table II contains 

the parameters of harmonic and biharmonic functions which have been htted to the data for the soliton marge and 
amplitude for three cases of the parameter K. Comparing CC theory with simulations, we see that the results for 
the mean values and the amplitudes of the first harmonics agree qualitatively. Finally we discuss the results for the 
translational motion of the solitons. There are very small oscillations of the soliton position q{t) around a mean 
trajectory v t, see Figs. |^and 10 We compute the discrete Fourier transform (DFT) of q{t) - vt for the CC theory 
and simulations and observe the same frequencies as above for the soliton charge and amplitude. 

Table H] shows that v is always close to the initial value vq = 0.1 which means that the translational motion is only 
weakly affected by the intrinsic soliton oscillations. The agreement between Vee and Vsim is not so good (the maximal 
error in Table |l] is about 14%). The reason is that the plane wave phonons with k = —K are not taken into account in 
the CC theory. The monotonic behavior of Vsim as a function of K is explained qualitatively in the following way: For 
positive K the phonon phase velocity is negative, which results in head-on collisions with the soliton. Here Vsim < vq, 
which can be explained by assuming negative spatial shifts of the soliton due to the collisions. For negative K the 
plane wave phonon overtakes the soliton which results in positive shifts explaining that Vsim > fo- 


VIII. SUMMARY 

We investigated how a solitary wave solution of the nonlinear Dirac (NLD) equation evolves in time under an 
external force with the components fj = rj eyip[—i{vjt — Kjx)]^ j = 1, 2. As an ansatz for a collective coordinate (CC) 
theory we took the exact Lorentz boosted solitary wave solution of the unperturbed NLD equation. The collective 
variables are the soliton position q{t), inverse width f3{t) and phase The variable /3 is related to the frequency 
uj{t) = a/ that appears in the solitary wave solution and lies in the range 0 < uj < m. In the non-relativistic 
regime ui is close to the mass m. The forced NLD equation is obtained in a standard way from a Lagrangian density. 
We restricted ourselves to the case with — 0, inserted our ansatz, integrated over space, and obtained the Lagrangian 
as a function of the collective coordinates. The Lagrange equations are three coupled ODEs. In two special cases we 
obtained approximate analytical solutions, but in general the ODEs were solved numerically by a MATHEMATICA 
program. We chose parameters and initial conditions from the non-relativistic regime where solitary wave solutions 
are expected to be stable. The solutions are periodic in time, which means that the solitary waves exhibit intrinsic 
oscillations with frequency flee- The translational motion is also affected, though weakly, because the position q{t) 












spectrum of Q{t) 
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FIG. 5. Intrinsic soliton oscillations for the case of constant, homogeneous force {K = 0) and zero initial velocity {vo = 0). 
Other parameters and initial condition (IC): g = 1, m = 1, r = 0.01, ojq = 0.9, (j)o = 7r/2, integration time tf = 800. Panels 
(a) and (b): charge from CC theory and simulation, respectively. Panels (c) and (d): amplitude a = 2[m — uj{t)]/g'^ and 
ma,Xx p{x,t), respectively. Panel (e): discrete Fourier transform (DFT) of Q{t), soliton peak at oji = 0.9032 and phonon peak 
at L 02 ~ 0.997 « Vl + with k = —K — 0. Panel (f): DFT of max^, p{x,t), peaks at loi = 0.9032, 013 = uj 2 — = 0.0942. 
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A'(7r/100) 

^CC 

^sim 

'^CC 

sim 

-6 

0.911062 

0.91106 

0.0855975 

0.10012 

-5 

0.911062 

0.91106 

0.0865441 

0.10008 

-4 

0.911062 

0.91106 

0.0914469 

0.10005 

-3 

0.903208 

0.90321 

0.101902 

0.10005 

-2 

0.903208 

0.90321 

0.111025 

0.10004 

-1 

0.895354 

0.89535 

0.114524 

0.10001 

0 

0.8875 

0.89535 

0.114569 

0.099968 

1 

0.8875 

0.89535 

0.112703 

0.099932 

2 

0.879646 

0.8875 

0.109807 

0.09987 

3 

0.879646 

0.8875 

0.10646 

0.099805 

4 

0.879646 

0.87965 

0.103058 

0.099745 

5 

0.871792 

0.87965 

0.0999051 

0.099677 

6 

0.871792 

0.87965 

0.0972293 

0.099607 


TABLE I. Frequency of intrinsic soliton oscillations and average soliton velocity as a function of the parameter K. We compare 
the results of the CC method and numerical simulation at various values of K in multiples of rr/lOO . Here we choose the initial 
conditions vq = 0.1, oiq = 0.9, (;Ao = -e I2,r\ = r = 0.01, r2 = A'2 = 0,tfin = 800, and denotes the integration time. 


Parameters 

Theory 

Simulations 

K = Q, vq = Q 

Q{t) « 0.97 - 0.072 cos(a;it -h 1.15) 
a{t) « 0.2 — 0.025 cos( 6 L)it -|- 1.15) 

Q{t) « 1.0097 — 0.049 cos(tuit -I- 0.058) — 0.034 cos(a;2t + 0.88) 
maxa; p « 0.2 — 0.0067 cos(cuit -|- 0.036) -I- 0.0065 cos(a;3t -|- 0.98) 

K = -37r/100, Vo = 0.1 

Q{t) « 0.97 - 0.073 cos(a;it -h 1.67) 
a{t) « 0.2 — 0.026 cos( 6 L)it -|- 1.67) 

Q{t) « 1.0092 - 0.071 cos(tJit-h 2.0239) - 0.038 cos(a;2t-6.74) 
maxa; p « 0.2 — 0.0096 cos(cuit -|- 2.02) -|- 0.0061 cos(tJ3t -|- 0.797) 

K = 37r/100, Vo = 0.1 

Q{t) « 0.97 - 0.0594 cos(u;it -f 1.31) 
a{t) « 0.2 — 0.021 cos(u;it -|- 1.31) 

Q{t) « 1.0093 — 0.054cos(a;it -I- 0.72) — 0.029 cos(a;2t — 0.2) 
maxj; p « 0.2 — 0.0089 cos(cuit + 0.72) — 0.00797 cos(a;3t — 0.26) 


TABLE 11. Least-squares fits to Q{t), amplitude a = 2[m — uj{t)]/g'^, and maxa, p from theory and simulation. Other parameters 
as in Figs. Note that 1 ^ 1 , uj 2 and wa are given in the captions of Figs. 



FIG. 6. Snapshots of the soliton profile at different times. Same parameters and initial conditions as in Fig. 
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(d) 




frequency frequency 

(e) (f) 


FIG. 7. Intrinsic soliton oscillations for the case of a harmonic inhomogeneous force {K = Svr/lOO ~ 0.094248) and vo = 0.1. 
Other parameters and initial conditions as in Fig. [q Panels (a) and (b): charge from CC theory and simulation, respectively. 
Panels (c) and (d): amplitude a = 2[m — Lij(t)\/ ^ axid maxi,p(x,t), respectively. Panel (e): DFT of Q(t), soliton peak at 
ui-i = 0.8875 and phonon peak at uj 2 = 1.0053 « \/l + with k = —K. Panel (f): DFT of max^, p{x,t), peaks at cui = 0.8875, 
a)3 = a;2 — car = 0.11781. 
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frequency frequency 

(e) (f) 


FIG. 8. Oscillations of the translational motion of the soliton. K = Stt/IOO « 0.094248 and vo = 0.1. Other parameters 
and initial conditions as in Fig. Panels (a) and (b): q(t) from CC theory and simulation, respectively. Panels (c) and (d): 
q(t) — vt from theory (v = Vcc = 0.1065) and simulation (v = Vaim = 0.099805), respectively. Panel (e): DFT of q{t) — Vcct, 
peak at cui = 0.87965. Panel (f): DFT of q{t) — Vsimt, peaks at uj\ = 0.8875, 2aJi = 1.7750 and at wa = a ;2 — = 0.11781. 
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(b) 





frequency frequency 

(e) (f) 


FIG. 9. Intrinsic soliton oscillations for the case of a harmonic inhomogeneous force {K = — Stt/IOO ~ —0.094248) and uo = 0.1. 
Other parameters and initial conditions as in Fig. [q Panels (a) and (b): charge from CC theory and simulation, respectively. 
Panels (c) and (d): amplitude a = 2[m — and ma.Xx p{x,t), respectively. Panel (e): DFT of Q{t), soliton peak at 

ui-i = 0.9032 and phonon peak at uj 2 = 1.0053 « \/l + with k = —K. Panel (f): DFT of max^, p{x,t), peaks at cui = 0.9032, 

U!3 = UJ2 — UJl = 0 . 1021 . 
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FIG. 10. Oscillations of the translational motion of the soliton. K = —Stt/IOO ~ —0.094248 and vq = 0.1. Other parameters 
and initial conditions as in Fig. Panels (a) and (b): q{t) from CC theory and simulation, respectively. Panels (c) and (d): 
q{t) — vt from theory (v = Fee = 0.1019) and simulation (v = Vsim = 0.10005), respectively. Panel (e): DFT of q{t) — Vect, 
soliton peaks at lo\ — 0.9032 and 2cni = 1.8064. Panel (f): DFT of q{t) — Vaimt, peaks at cJi = 0.9032 and 1.7907 « 2a;i. 
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oscillates around a mean trajectory Vcct- We compared our CC predictions with numerical simulations of the forced 
NLD equation: The solitary wave solutions are in fact stable and periodic. The observed frequency ^sim = is 
nearly identical with flee- However, Vsim agrees with v^c with an error of about 14%. The reason for this is that the 
CC theory does not include phonons (short for linear excitations). In fact, a specific plane wave phonon mode with 
wavenumber k = —K is excited together with the intrinsic oscillations in order to conserve the total momentum. The 
predicted frequency fix = 'Jrn? + agrees perfectly with the frequency u }2 in the spectra of all variables. 

For the future work we plan to take initial conditions away from the non-relativistic regime, i.e. initial uj not close 
to TO, and initial velocity not much smaller than the speed of light. Moreover, it will be very interesting to see what 
is the influence of time dependent external forces, i.e. non-vanishing Vj. 
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Appendix A: Relevant Integrals 


For our ansatz in the rest frame, we have that for k = 1 the charge Q is 


Q= f = f dxiA^ + B^) = . f 


, I J- tanh^ fdx , , ^ 2/3 

dx~ ---sech px = 

" 5 w 


g‘^{m + uj) J_ao (1 — q:^ tanh^/3a:)2 
For Sec. V we need explicit expressions for the following integrals (in what follows, y = tanh/3a;): 


Hi= dx - 9,^4-7^4'] = / dx{B'A - A'B) = 


2 (to — (jj) 


a dy 


l-J/2 


-1 (1 - a^y'^Y 


2 tanh 


m — u! 
to + m 


-0]= In 


(Al) 


(A2) 


H 2 = m J = mil = y dx{A^ — B^) = 

4to , _i , , 

= tanh ^(a) = Mq, 

r 


2to/3 


dy 


4to/ 3 tanh ^(a) 


g‘^{m + ui)J_i (1 — Qf^y^) g‘^{m + uj) 


where Mq is the mass in the rest frame. Note that Mq has the property of going to zero as w —)■ 1. 


I2= j dx{A^ - B^Y = ^ f 

J g^{m + ujY J- 


dy 


I-z/2 


4^3 


(a^ + 1) tanh ^(a) — a 


g'^{m + ujYJ-i {1 — a'^y^Y g'^im + uiY 


2 2 

= 4^0 = 4 ^ 1 - 

g g 


To calculate the integral Jj defined in (5.181, first we rewrite it as 


r+oo 


Jj{u:,q) = / dzA{z) cos{20ajz) = 


yj2{m + w)/3 


^ + 00 


dz 


guj 


cosh(/3z) cosh{i20ajz) 
— -I- cosh(2/3z) 




TO -I- u:)0 cosh[(l J- i2aj)l3z] + cosh[(l — i2aj)0z] 

— / dz - 


gw 


■ cos\i{20z) 


(A3) 


(A4) 


(A5) 
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Now using expression (6) on page 357 of [5Tj, after some straightforward calculations we obtain 


Jj(w,g) = 


TT COS bj 


gy/uj coshujir’ 


where Uj and bj are defined in Eq. (5.181. The integral Nj can be calculated in a similar way. 
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